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We visualize the Kohn-Sham kinetic energy density (KED), and the ingredients - the electron density, its 
gradient and Laplacian - used to construct orbital-free models of it, for the AE6 test set of molecules. These 
are compared to related quantities used in metaGGA’s, to characterize two important limits - the gradient 
expansion and the localized-electron limit typified by the covalent bond. We find the second-order gradient 
expansion of the KED to be a surprisingly successful predictor of the exact KED, particularly at low densities 
where this approximation fails for exchange. This contradicts the conjointness conjecture that the optimal 
enhancement factors for orbital-free kinetic and exchange energy functionals are closely similar in form. In 
addition we find significant problems with a recent metaGGA-level orbital-free KED, especially for regions 
of strong electron localization. We define an orbital-free description of electron localization and a revised 
metaGGA that improves upon atomization energies significantly. 
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I. INTRODUCTION 

The Kohn-Sham kinetic energy density (KED) - the 
kinetic energy per volume defined by the orbitals gener¬ 
ated by the Kohn-Sham equation - plays a central role 
in the development of density functional theory (DFT). 

In the “Jacob’s Ladder” paradigm for characterizing the 
exchange-correlation (XG) energy in density functional 
theory,Iil the KED is the key variable of the central, 
metaGGA rung of functionals.l^^ As a local energy den¬ 
sity, it provides information about electronic structure 
complementary to that provided by the local electron 
density and its derivatives that describe lower rungs of 
DFT. Particularly important is its ability to distinguish 
between regions of electron localization,^"!! for which self¬ 
interaction error is important, and regions of delocaliza¬ 
tion such as metals where they are not. 

The centrality of the KED in DFT development is 
highlighted by the implicit role it plays in the rungs of 
DFT lower than the metaGGA. These may be thought 
of as a Jacob’s ladder of approximations to the KED 
as much as one of approximations to the exchange- 
correlation energy. The lowest rung of the XG ladder, 
the local density approximation orLDA,!^ corresponds to 
the Thomas-Fermi approximatiorP^*^ to the KED. The 
more commonly used generalized gradient approximation 
(GGA) introduces, in addition to the local density, the 
gradient of the density as a variable in XG functional 
construction. The same information is contained in the 
von Weizsacker KEEC^ that describes the KE of local¬ 
ized electron pairs. It has been used to generate a large 
number of G GA’s for the KED, both empiricaP^l^ and 
nonempirical,!J^^*J^ though not with the success they have 
enjoyed in describing the XC energy. To describe the XG 
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energy at the next, metaGGA, level of the theory, not 
only the KED, but also the Laplacian of the densitjff^lf^ 
may be used as an additional variable in functional con¬ 
struction. A metaGGA descri ption o f the KED is thus 
possible, using the Laplacian.^SESEl similarity of 

the KED and XC functional ladders leads to a conjoint¬ 
ness conjectur^^that the optimal orbital-free correction 
to the Thomas-Fermi KE is similar in form to that for 
LDA exchange. More importantly, it enables one to ap¬ 
ply lessons learned in constructing the one functional to 
constructing the other. This is important for orbital-free 
DFT, in which the Kohn-Sham KED is replaced by an 
explicit functional of the density, removing completely 
the need for orbitals. 

The orbital-free modeling of the K ED h as taken on in¬ 
creasing importance in recent years.^^fiHli Given a cubic 
scaling in the number orbitals, the Kohn-Sham method 
becomes prohibitively expensive for large-scale appli¬ 
cations that require the accuracy of atomistic simula¬ 
tion. These involve applications such as the dynam¬ 
ics of nanoscale material^^ requiring intrinsically large 
system size, high throughput as in alloy design, or a 
need for large number of excited states, as can occur for 
finite-temp eratur e applications such as warm dense mat¬ 
ter (WDM) .[sum A completely orbital-free density func¬ 
tional theory (OFDFT), using an orbital-free expression 
for the kinetic energy, becomes an important tool in these 
cases. Unfortunately, OFDFT is inherently less accurate 
than Kohn-Sham DFT; for example, the Thomas-Fermi 
approximation is unable to predict molecular binding,!^ 
something the LDA has no problem doing. Nevertheless, 
nonlocal model^^^^^^ have achieved reasonably high ac¬ 
curacy, allo wing f or impressive calculations for solid-state 
applicationsalbeit within the limitation of requiring 
different functionals for different material classes. And, 
with a focus on improving potentials and thus forces in 
the context of WDM, a number of GGA-level OFDFT’s 















have been developed in recent These co¬ 

incide with impro vements in infrastructure for practical 
calculations P^ I ^^ ^ 

A third role of the KED has been as a tool for vi¬ 
sualizing the electronic structure of the chemical bond. 
The kinet ic ene rgy density has been the subject of in¬ 
vest igatioiP^Eni particularly as a localized-orbital locator 
(LOL)j^i^and an impressive number of related quanti¬ 
ties have been defined and investigated as wellPerhaps 
the most popular is the electron localization factor ELF, 
which is based upon a comparison of the Kohn-S ham to 
Thomas-Fermi and von Weizsacker KED’s-^^^^^^^ It is of 
particular importance for the development of metaGGA’s 
for the XC energ y an d in the conceptual understanding 
of why they work.l^ Also of note is the quantum theory 
of atoms in molecules (QTAIM),I 2 SH 1 II an approach to vi¬ 
sualization which is in some ways an orbital-free version 
of ELF analysis, using gradients and Laplacians of the 
density to analyze bonding structures. 

Despite the strong connecti on betw een the arguments 
used to build XC metaGGA’P®^!^ and those used to 
visualize the chemical bond, the tool of visualization has 
not often been used to provide feedback into DFT de¬ 
velopment. The properties of the exchange-correlation 
hole, describing the hole around an electron caused by 
Pauli exclusion and Coulomb repulsion have been an 
important tool in the construction of both GGA’s and 
their successors.^^^^^ In particular, the visualization of 
the hole has be en a valuable tool in assessing the ac¬ 
curacy of DFT’s. 1 ^ 2 *^ However, the exchange-correlation 
hole is a difficult many-body calculation, and the depen¬ 
dence of measurables like the atomization energy or bond 
length on the nature of the XC hole occurs implicitly 
through the mediation of complex functionals and thus 
is hard t o det ermine. (But connections can sometimes 
be made.^^^^^ In the case of the KED, however, visual¬ 
ization is of direct helpPMl 

- how an orbital-free density 
functional theory for the Kohn-Sham KED actually com¬ 
pares to the real thing requires no more than running a 
standard DFT code and visualizing the results. 

In this paper, we perform highly converged Kohn- 
Sham DFT calculations and visualize the electron den¬ 
sity, its gradient and Laplacian, the KED and some ap¬ 
proximations for these used in DFT, for the AEG test set 
of molecules, in a pseudopotential plane wave approach. 

The AEG test sefP^ is a set of G molecules - Cyclobu¬ 
tane (C4H8), Propyne (C3H4), Glyoxal (C2H2O2), Sili¬ 
con Monoxide (SiO), Disulfur (S2), and Silicon Tetrahy- 
dride (SiH4) - chosen for their ability to reproduce the 
average atomization energy of common DFT’s over much 
larger test sets. For such a small set the AEG shows a 
richness of bonding scenarios - single, double, and triple 
bonds, covalent to nearly ionic, including first and second 
row atoms, and a large-cation, small-anion system simi¬ 
lar to important semiconductors like GaN. Thus it covers 
many situations commonly seen in organic chemistry and 
in semiconductors as well. 

Our motive for using pseudopotentials is two-fold. 


First of all, many current O FDFT applications rely upon 
the use of pseudopotentials,l^2*^ although more accurate 
approaches do exist . 1 ^ More importantly, the pseudopo¬ 
tential plane-wave approach permits an arbitrary con¬ 
vergence of the particle density associated with the pseu¬ 
dopotential and thus a map between a u—representable 
density and the related KE density that is as accurate 
as possible. It thus gives insight into the universal map 
between kinetic energy and density that is a corollary 
of the Hohenberg-Kohn theorem. Although the method 
does not produce the correct density for real molecules, 
and thus introduces errors into the chemical characteriza¬ 
tion of the test set, it arguably gives us simpler problem 
to model, and much of what is learned for pseudopoten¬ 
tial systems should help to construct functionals for the 
all-electron case.^^ The use of pseudopotentials enables 
particularly the study of asymptotic features not possible 
with a typical gaussian basis set. 

Finally, the choice of exchange-correlation functional is 
irrelevant to the universal mapping between the Kohn- 
Sham kinetic energy and the charge density, in which 
the electrostatic potential energy plays no role. We work 
with the LDA and PBE exchange-correlation functionals, 
which produce reasonably accurate bond lengths for the 
test set and should produce densities and orbitals close 
to the exact ones for pseudopotential systems. 

In our visualization, we have deemphasized (but do not 
ignore) the ELF, already studied extensively for a large 
number of molecular systems. We look rather at the basic 
ingredients of the orbital-free KED, the electron density 
n, and related derivatives |Vn|^ /n and V^n, focusing es¬ 
pecially on applications of their use in DFT. One is a 
common approximation based upon the gradient expan¬ 
sion in the limit of slowly varying densities used in many 
metaGGA’s to replace V^n, a natural descriptor in this 
limit, for tks- The second is a sophisticated metaGGA- 
level orbital-free model of the Kohn-Sham KED, the 
mGGA.l^ This takes advantage of lessons learned in de¬ 
veloping metaGGA’s for exchange, particularly of defin¬ 
ing and respecting key constraints and limiting cases for 
the kinetic energy. Despite the promise of its design phi¬ 
losophy, the mGGA has deficiencies - its potential does 
not bind moleculeJi^ and even used non-self-consistently 
fails to improve upon Thomas-Fermi predictions of at¬ 
omization energies .^However, it is of value as a starting 
point of thinking how to construct a metaGGA; and since 
it is a model of the kinetic energy density as such, it is 
directly testable by visualization of this quantity. Our 
investigation of the mGGA shows, despite its excellent 
description of atomic KED’s, surprising failures in its de¬ 
scription of the KED of bonds, and thus in its prediction 
of atomization energies. Our visualization work makes it 
easy to diagnose and suggest a fix to this problem, one 
which defines, and demonstrates at least in an ad hoc 
fashion, a potential lower bound to the KED. 

The rest of this paper is organized as follows: Sec. |n] 
describes the theoretical background of the paper - the 
density functional theory of the kinetic energy and its 
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relation to exchange in metaGGA’s. Sec. Ill covers the 
basic methodology used. Sec. |T^ details the chief results 
of visualization, while Sec. |V] applies the lessons learned 
to construct and make preliminary tests of a correction 
to the Perdew-Gonstantin mGGA and Sec. VI presents 
our conclusions. 


This is the tak e-off point for many orbital-free function¬ 
als for point of view consid¬ 


ered in this p aper. At the same time, nonlocal function- 
aipDEMIIll take the form 


TksVA= J j n°‘{v)W{v,v')n^{r')cPrcf'r' 


( 6 ) 


II. THEORY 

The positive definite form of the kinetic energy density 
in Kohn-Sham theory is given by 


which may be related to the semilocal picture through 
an expansion of the kernel W(r, r') for small r — r'.l^ 
The lowest level of semilocal functional - the equivalent 
to the LDA in XG functionals - is the Thomas-Fermi 
model. 


occup 

Tics= 2 E ( 1 ) 

i 

where (j)i are Kohn-Sham orbitals from which the electron 
density is constructed: 


occup 

= E 

i 

and fi is the occupation number of each orbital. Integra¬ 
tion over all space gives the kinetic energy 

TifsN = J TKs{T^)d^r. ( 3 ) 

A generalization in terms of the spin density and spin- 
decomposed KED’s is easily constructed by restricting 
the sums in the equations above to a specific spin species 
but will not be considered here. The KED is well defined 
only up to the arbitrary addition of a divergence of a 
vector function - the integration of such an addition is 
zero and leaves the physical measurable Tks unchanged. 
Thus any number of physically equivalent KED’s may be 
constructed, with a common alternative to Eq. Q being 

occup ^ 

^'ks = -^Y. ( 4 ) 

i 

The value of Eq. Q is that it is positive-definite like the 
particle density, and that a number of properties of the 
KE are conveniently framed in terms of it. 

The key principle for this paper is that since Tks [n] is 
a functional of the ground state electron density n, tks 
must be one too. There exists some map TKs[n] from 
Eq.© to Eq. (|lj) that need not explicitly rely on orbitals. 
However the form of this map is unknown, and unlike the 
exchange-correlation functional of standard Kohn-Sham 
theory, approximate functionals are often far from satis¬ 
factory. Specifically, as is done in the lower rungs of the 
XC ladder of approximations, we can define a “semilo¬ 
cal” model of TKs['n]^ in terms of functions of the local 
density and its derivatives: 

T7J’™"[n]= y’T“PP™"[n(r),Vn(r),V2n(r)]d3r (5) 


ttf = ( 7 ) 

with kp = (STT^n)^/^ the fermi wavevector of the ho¬ 
mogeneous electron gas (HEG). At the next level of 
appro ximat ion, the gradient expansion approximation 
(GEA)£S 16 ZI of the KED is given by: 

tgea = ttf + + 0{W^). ( 8 ) 

Terms up to fourtlP^ and sixth ordeil^ in this expansion 
are known. 

As is the case with exchange, in order to preserve the 
proper scaling of Tks under the uniform scaling of the 
charge density, the form of an orbital-free functional for 
the KED is restricted to that of a function of scale- 
invariant quantities times the local density approxima¬ 
tion. Thus the GEA can be recast as 


Tgea = 




Ttf, 


( 9 ) 


in terms of invariant quantities: 


|Vn|^ 

^ AkpU ’ 

(10) 


(11) 

^ Aklpn 


Similarly, the most general form for a semilocal func¬ 
tional is a generalization of the GEA form in terms of an 
enhancement factor Fg modifying ttf'- 


'Tsemilocal = Fs{p,q)TTF- ( 12 ) 

The enhancement factor Fs for the kinetic energy plays 
a similar role to that for exchange, Fx, in conventional 
GGA’s, where the exchange energy density is expressed 
as a correction to the LDA in the form FxCx^^- 

In constructing generalized gradient functionals, it is 
conventional to omit the term proportional to in 
the GEA expansion as this integrates to zero and does 
not contribute to the overall kinetic energy.!^ Then, by 
approximating the gradient expansion to all orders in the 
remaining variable p, one obtains the next natural step, 
the GGA.Ii^ However, our goal is to visualize the local 
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quantity rifs(r), and for this purpose, the term in its 
gradient expansion cannot be ignored. Moreover, keeping 
it is necessary to implement local constraints on orbital- 
free approximation to tks (and thus constraints on Tks) 
correctly, and we do so in the work that follows. V^n 
is normally considered as a higher-order variable whose 
introduction in a functional defines the next, metaGGA, 
rung of functionals. 

Up to this level of approximation, the process of build¬ 
ing a kinetic energy functional mirrors that for exchange, 
so that the conjointness conjecture has been mad^^that 
the optimal form for each functional at a given level of 
approximation are closely related: Fs ~ Fx- This re¬ 
lationship has never been explicitly defined, but is nor¬ 
mally taken to be that of nearly identical functional forms 
with different constants.l^^®^ This strict conjecture has 
been demonstrated to be wrong,!^ but a philosophy of 
conjointness, using the experience of designing exchange- 
correlation functionals to i nform th e design of KE func¬ 
tionals, is common practicel ^^ l ^^ l ^'^ l 

Nevertheless, there are fundamental differences be¬ 
tween the two functionals, particularly in the physics 
of the large inhomogeneity limit p,\q\ S> 1. For the 
Kohn-Sham KED in real systems, the most crucial is¬ 
sue is the limit of a one-particle system or two-particle 
spin-singlet system. In this case it reduces to the the 
von Weizsackeiti^ functional: 


TvW = o 


1 |Vn|' 


8 n 


(13) 


This is the exact result for a system of N particles obey¬ 
ing Bose statistics, so that in the ground state they oc¬ 
cupy a single ground state orbital, (j)o = s/n/N. The KED 
needed to create the density n(r) with fermions, that is, 
the energetic cost of Pauli exclusion, is given by the dif¬ 
ference between the Kohn-Sham and Wigner KED’s and 
must be positive definite!^ 


TPauli = Tks — TvW > 0. 


(14) 


Notably, this von Weizsacker lower bound is not re¬ 
spected by the GEA. If we rewrite Eq. (13) in terms of 
an enhancement factor, we find Fg^ = 5p/3 - a depen¬ 
dence on p that is nine times faster than that of the GEA. 
For (7 = 0, TOE A falls below Tvw for the relatively mod¬ 
est value of p = 27/40. The constraint can be imposed 
by changing the coefficient in the gradient expansion to 
5/3, in which case the slowly-varying limit is incorrect. 
In contrast, exchange is constrained by the Lieb-Oxford 
bouncP^ that limits the contribution from the low den¬ 
sity tail outside the classically allowed range of electron. 
This limit has no intrinsic tie to the single-orbital limit 
and we shall see that the KED behaves very differently 
from exchange in this limit. 

Recently a metaGGA-level KED functional of the form 
of Eq. ( [l^ , the Perdew-Constantin rnGGA,!^ has been 
developed by applying lessons learned in constructing 
constraint-based exchange-correlation functionals. It sat¬ 
isfies the gradient expansion up to fourth-order in the 


limit of slowly varying density and the von Weizsacker 
bound and other constraints for large values of p and 
q. The function interpolates between the gradient ex¬ 
pansion and von Weizsacker limits using a nonanalytic 
but smooth interpolating function that depends on an 
effective localization measure z = Fqea-m — Fyw , with 
Fqea-m a metaGGA designed to be the best-possible 
analytic functional built from the starting point of the 
slowly varying electron gas. This is explicitly a model of 
the kinetic energy density, designed to take the place of 
the KED in exchange-correlation metaGGA’s, and thus 
is meaningfully tested by means of visualization of the 
KED. 

So far, the Jacob’s Ladder of approximations of the 
Kohn-Sham KED parallels the development of exchange 
functionals. A divergence now occurs in that, for ex¬ 
change and correlation, the KED itself can be used as 
a variable for building further approximations. In the 
standard approactP to constructing metaGGA’s for ex¬ 
change, the Laplacian of the density, which appears ir- 
reducibly in fourth and higher-order terms in the gradi¬ 
ent expansion is introduced implicitly through the use of 
the Kohn-Sham KED. This is achieved by rewriting the 
gradient expansion for tkSi [Eqs. (|8j|^], to construct a 
replacement for V^n, good to second order in this expan¬ 
sion. This “pseudo-Laplacian” is given by: 

= 6(rifs - ttp’) - ^|Vnp/n, (15) 

which then replaces V^n in the construction of the 
metaGGA. V^n approaches in the limit of slowly 
varying density, deviating from it only where the V^n 
gets large, such as at the cusp in the electron density at 
the nucleus. It is unknown how well this approximation 
works in practice for features of electronic structure like 
covalent bonds, which locally may have small p and q but 
are part of systems that are far from the slowly-varying 
limit globally. This quantity can then serve to test the 
quality of the GEA. 

Perhaps the most physically significant role played by 
the KED in a metaGGA is as a measure of electron 
localization.!^ This is done by taking the ratio of the 
Pauli contribution to the Kohn-Sham KED to that of 
the Thomas-Fermi model. 


a = 


Tks -Tvw 
ttf 


(16) 


In regions where the KE density is determined predomi¬ 
nantly by a single molecular orbital, tks approaches Tvw 
and a —)■ 0. This limit describes single covalent bonds 
and lone pairs, and generally situations in which the self¬ 
interaction errors in the GGA and LDA are most acute. 
The HEG, and presumably systems formed by metallic 
bonds, corresponds to tks = ttFj Tvw 0 and a ^ 1. 
Between atomic shells and at low density, 1, poten¬ 
tially tending to oo for an exponentially decaying density 
if Tpauli vanishes more slowly than This limit can 
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be used to detect weak bonds such as van der Waals in¬ 
teractions and define interstitial regions in semiconductor 
systems. The information on local environment can then 
be used to customize gradient approximations for specific 
subsystems.^ 

It is a short step from a to the electron localization fac¬ 
tor or ELFI^ used in the visualization of electronic struc¬ 
ture: 

ELF=-^. (17) 

1 - 1 - 

This converts the information contained in a to a func¬ 
tion between one, when a = 0, to zero (a oo), useful 
for visualization, but less so in functional construction. 
The related LOlP is closer in form to V^n, and is basi¬ 
cally the enhancement factor Fks='^ks/ttf recast into 
a convenient form: LOL = 1/(1 -|- F^s)- 

Finally we note that the a used in defining the ELF is 
also the enhancement factor for the Pauli KE: TpauU = 
aTTF- And thus one can consider the project of con¬ 
structing OFDFT as intimately tied to the project of vi¬ 
sualizing electronic structure - constructing orbital-free 
models to the ELF and the information on electron lo¬ 
calization it contains. This has b een t he perspective of 
several recent studies of the 


III. METHODOLOGY 

As noted in the introduction, we use the plane-wave 
pseudopotential method for performing DFT calculations 
- this allows us to solve nearly exactly the Kohn-Sham 
equation for a model system and acquire highly accurate 
orbitals, but for an approximate system. For this pur¬ 
pose, the ABINIT plane-wave pseudopotential code^SE^ 
was employed with an LDA and PBE XC functionals. 
Standard Troullier-Martins pseudopotential^^ from the 
ABINIT library were used for both. Geometries were op¬ 
timized using the Broyden-Fletcher-Goldfarb-Shanno al¬ 
gorithm, ESI to a force tolerance of 5 x 10“^ hartree/bohr. 

The main convergence error in our calculations was 
that of using a finite-sized periodic simulation cell, ne¬ 
cessitated by the use of a plane-wave expansion. The 
simulation cell size was chosen so that total energies were 
converged to within 3 x 10“® hartree. Errors in nearest- 
neighbor bond-lengths due to finite system size are less 
than 5 x 10“® A. In order to get good spatial resolu¬ 
tion of plots, we took a plane-wave cutoff of 99 hartree 
for all systems, well above that needed for convergence 
of energies to chemical accuracy (< 40 hartree) in the 
pseudopotential systems. The convergence errors in to¬ 
tal energy from the finite plane-wave cutoff range from 
10 “^ hartree for SiH 4 to I 0 “® for C 4 H 8 and the error 
in nearest-neighbor bond-lengths, from 10“^ to 10“® A. 
Converged simulation-cell parameters for each molecule 
may be found in the supplementary information. 

Given a periodic cell, the density and related expec¬ 
tations should suffer boundary effects. Most notably. 


whereas the density and its derivatives and the kinetic 
energy density should decay exponentially to zero in a fi¬ 
nite system, these will approach a small finite value at the 
cell boundary. For the cell sizes used, this minimum value 
of the density is on the order of 10“® a.u., for systems 
with maximum densities on the order of an a.u.; a sig¬ 
nificant distortion from exponential decay is observable 
only within two bohr of the location of the minimum. 

The ABINIT code outputs density and KE density as a 
three-dimensional uniform grid over the periodic simula¬ 
tion cell, with grid spacing determined by the dimensions 
of the fast Fourier transform used in the plane-wave code. 
The real-space grid used to accommodate a 99 hartree 
plane-wave cutoff has a resolution of 0.11 bohr, defining 
the resolution of our plots. The Laplacian and gradient 
of density were evaluated numerically on this grid using 
a Lagrange-interpolating polynomial method. Color sur¬ 
face plots and contours were generated using gnuplot and 
the associated pm3d utility. 

IV. RESULTS 

First of all, to assess the quality of data within the 
plane-wave pseudopotential approach, we show results 
for basic structural properties for the AE6 test set. Ta¬ 
ble |l] shows the mean relative error (MRE) and mean 
absolute relative error (MARE) of LDA and PBE pseu¬ 
dopotential predictions of bond lengths for the AE6 test 
set, as compared to experimental data. The LDA gives 
an excellent fit to double and triple bonds and about a 
1% over-binding of singl e bonds, in line with other re¬ 
sults for the LDA .1^212218^ An atypical tendency to under¬ 
bind for C-H bonds leads to a MRE whose accuracy we 
suspect would not hold for larger test sets. The overall 
tendency of the PBE is to increase bond lengths relative 
to the LDA, again the expected trend, which results in a 
slightly better absolute agreement with experiment. 



LDA 

PBE 

MRE (%) 

-0.006 

-0.087 

MARE (%) 

0.68 

0.53 


TABLE I. Performance of pseudopotential DFT calculations 
for the bond lengths of the AE6 test set - mean relative er¬ 
ror (MRE) and absolute relative error (MARE) in angstroms 
compared to experimental data from Ref. M 

The summary performance of DFT predictions for at¬ 
omization energies is shown in Table [ll| Again the trend 
of the LDA is to over-bind with respect to experiment 
and that of the PBE to remove much of this error. The 
LDA does worst energetically for systems with a double 
or triple bond: S 2 , SiO and C 2 H 2 O 2 . Our pseudopo¬ 
tential estimate of the MAE for the PBE functional on 
the AE6 test set compares reasonabl y wel l with those ob¬ 
tained from all-electron calculationa^^ISSl using gaussian 
basis sets. The two approaches agree for singly-bonded 
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LDA PBE PBE-ae 

MSE 

67.1 

20.8 

12.0 

MAE 

67.1 

23.0 

15.1 

MARE (%) 

16.0 

7.4 

4.4 


TABLE II. Summary errors (mean signed, mean absolute 
and mean absolute error in percent) of pseudopotential DFT 
calculations, and of an all-electron PBE calculatiorl^ for the 
atomization energy of the AE6 test set, measured relative to 
experimental data from Ref. isni In kcal/mol. 

systems while our pseudopotential approximation overes¬ 
timates the atomization energy of molecules with double 
bonds by about 10 kcal/mol per double-bond. A purely 
numerical calculation on an ultrafine gricP^ reports a 
MAE of 3.0 kcal/mol per bond for the PBE functional 
as compared to 3.6 kcal/mol per bond here, indicating 
that use of a pseudopotential overestimates binding but 
perhaps not by as much as indicated by the gaussian 
basis-set calculations. In any case, this error is minute in 
comparison to the large errors between orbital-free and 
Kohn-Sham kinetic energies. 

Further information about the convergence with re¬ 
spect to the finite size of the cell is shown in the supple¬ 
mentary material for this paper,including converged 
finite-size cell parameters for each molecule in Table S-I 
and finite cell boundary errors for S 2 in Fig. S-1. Ta¬ 
ble S-II shows per-molecule data from LDA and PBE 
pseudopotential calculations of the bond lengths of the 
AEG test set, compared to experimental data, and S-III 
does the same for atomization energies. 


A. Electronic structure: atoms 



r(a.u.) 


FIG. 1. (color online) Comparison of all-electron and 
pseudopotential kinetic energy densities for the carbon atom. 
Shown are the radial probability density versus distance from 
nucleus for the Thomas Fermi (solid line), von Weizsacker 
(dashed), and Kohn-Sham KED (lighter solid) as well as the 
quantity u^w (dot-dashed) defined in the text. The equiva¬ 
lent pseudopotential quantity for each is shown as a dotted 
line, matching at the cutoff radius 1.498 as- 

on is quite significant: the region of peak negative 
(the valence shell charge concentration or VSCC 
in QTAIM analysis) is broader in extent and the posi¬ 
tion of the critical point about 40% closer to the nucleus 
than in all-electron calculations. The maximum negative 
value of V^n is typically three times larger than its all¬ 
electron equivalent, with similar errors for the VSCC’s 
of molecules. Plots shown below for V^n and tks in 
molecules do faithfully reproduce the qualitative topo¬ 
logical features of the all-electron case, and are quantita¬ 
tively accurate at bond centers and asymptotically; but 
they must be treated with caution with respect to other 
quantitative details. 


Before showing results for molecules, it is instructive 
to compare pseudopotential and all-electron results for 
atoms. Fig. [^demonstrates this comparison for the den¬ 
sity, its gradient and Laplacian and the Kohn-Sham KED 
of the C atom. In order to make a clean comparison be¬ 
tween quantities, we convert the first three functions into 
equivalent kinetic-energy density models: ttf 
T yw = |Vn|^/ 8 n, and Uyw = The last is gen¬ 

erated by taking the functional derivative of Tyw with 
respect to density. 

The pseudopotentials for C are designed so that the 
pseudo-valence density matches the all-electron density 
after a cutoff radius of 1.498 gb- This match is respected 
for the other quantities as well. However the pseudoden¬ 
sity and thus ttf continue to match the real density with 
reasonable agreement almost to the core-valence transi¬ 
tion radius at about 0.8 as- The other pseudo-quantities 
deviate from their all-electron equivalents much more 
quickly, especially tks and Uyw- In the all-electron case, 
core orbitals smooth out the density and thus reduce the 
magnitude of the < 0 peak in the valence shell, and 
they add extra terms to txs- The quantitative impact 


B. Electronic Structure: molecules 


Figures 1^1^ 13 and [3 show contour plots for the den¬ 
sity and related quantities for pseudopotential models of 
several of the molecules of the test set: C 3 H 4 , C 2 H 2 O 2 
and SiH 4 . and SiO. In Fig. we show in the first row 
(a) the ground-state pseudo-density n and (b) the gra¬ 


dient factor |Vn|^ /n that appears in the gradient ex¬ 
pansion of Tks [Eq. (|3)] and the von-Weizsacker KED 
[Eq. (13l]. The second row shows (c) the Laplacian 
of the density V^n, and the gradient-expansion derived 
pseudoLaplacian [Eq. (15l] used in metaGGA’s. The 
third row shows (e) the Kohn-Sham KED tks and (f), 
the Perdew-Constantin mGGA model for the same. All 
quantities are plotted in hartree atomic units; all except 
(a) are thus dimensionally energy densities. The other 
figures show subsets of this suite of data, as identified by 
subcaptions, for the other three molecules. 

Each figure shows a two-dimensional slice through the 
molecule, with a color surface plot with values ranging 
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from blue (minimum value shown) to red (maximum). 
The numerical scales for the surface plots are shown in 
the bar to the right of each subplot. Superimposed upon 
these are contour plots. Thicker contour lines for the 
Laplacian and pseudo-Laplacian [Fig. |^c) and (d)] in¬ 
dicate the zero contour; the other four functions plot¬ 
ted are positive definite. The contours are adjusted to 
bring out details of bonding regions, and do not cover 
the atom cores. Contour values and ranges for the equiv¬ 
alent quantities of Laplacian and pseudo-Laplacian are 
identical, as are those for the two KED’s in subplots (e) 
and (f). Atoms and bonds in the plane of a plot are indi¬ 
cated by black dots and thick black lines; projections of 
out-of-plane atoms and bonds onto the plane of the plot 
are shown as open circles and thick dashed lines. 

We start out with a discussion of the four direct mea¬ 
sures of electronic structure - the density, | Vn|^ /n, V^n, 
and r and consider the two approximations V^n and 
TmGGA in following subsections. 

The subject of Fig. C 3 H 4 , is perhaps the richest ex¬ 
ample of the test set, illustrating several types of bonds. 
Three hydrogens (H^ to H^) bond with a tetrahedral ge¬ 
ometry to a carbon (C^), which is joined to the second 
carbon (C^) through a single bond, the second carbon 
shares a triple bond with the third (C^), and this is ter¬ 
minated with a final carbon-hydrogen bond. Our plots 
show a cut through the three central carbon pseudo¬ 
atoms aligned along the a;-axis and one hydrogen on ei¬ 
ther end; the other two hydrogens extend out on either 
side of the far left-hand side of the plane. 

The valence particle density (a) shows some features 
common to each molecule of the set. As in the single¬ 
atom case (Fig. [^, the density tends smoothly to a min¬ 
imum at the center of each carbon pseudo-atom. In con¬ 
trast, the hydrogen atom has no core electrons and the 
effect of the pseudopotential is simply to smooth out the 
cusp in the density at the nucleus. The highest elec¬ 
tron density is thus naturally within bonds - especially 
the triple bond (C^ — C^). The gradient-squared of the 
density (b) is nearly zero in regions with bonding, where 
the Laplacian (c) shows most structure, and is largest in 
the pseudo-atom core and at the edges of the molecule 
where is zero, as indicated by the thicker contour 
line in (c). The Laplacian is negative almost entirely 
along the center except for the interior of each carbon 
pseudo-atom. T his is a hallmark of covalent bonding in 
QTAIM analysiJ^^^^ - the center of a bond is a sad¬ 
dle critical point for the particle density, with a negative 
value for a covalent bond because of the the buildup of 
charge between atoms. The value of V^n at the C^-C^ 
bond critical point is -0.700 and that for C^-C^ is -1.143, 
reasonable values for C-C bonds.^^ V^n is positive in the 
pseudo-atom core, where the density is at a local mini¬ 
mum, and in the classically forbidden region far from the 
molecule. The Kohn-Sham kinetic energy density (e), is 
the smoothest and least structured of the measures of 
the Kohn-Sham system shown. As its relationship to the 
electron density is nontrivial, it not surprisingly appears 


to have little apparent correlation with it. It is primar¬ 
ily concentrated in the pseudo-atom cores with a strong 
peak at the center of the pseudo-atom. This follows the 
qualitative trend of the KED of all-electron systems,!^ 
except for the absence of shell structure. Otherwise it is 
significantly larger in the triple bond than in the single 
bonds, where it is nearly zero. 


Fig.0 shows Laplacian and KED quantities for the 
C 2 H 2 O 2 pseudo-molecule. This has a trigonal-planar 
form with a line of symmetry through the center of the 
molecule. The oxygens share polar double bonds with 
the carbons, and the carbons form covalent single sp^ 
bonds with each other and the hydrogens. There are two 
lone pairs of electrons present on each oxygen. The plot 
shows one oxygen, carbon and hydrogen, and part of the 
C-C bond at the bottom of the plot. The single C-C and 
C-H bonds are very similar to those in C 3 H 4 , so that the 
scale is adjusted to favor the oxygen atom which has a 
much larger density and KE density. 

The C-0 bond, being polar, exhibits several features 
not seen in C 3 H 4 . The density gradient is nonzero in the 
bond - the push of density towards the more electroneg¬ 
ative oxygen causes a local saddle point in the gradient 
on the oxygen side. VSCC lobes due to two sets of un¬ 
paired electrons are identifiable on the oxygen, but none 
on the bond axis, a reflection of the change in character 
of the bond. However, the VSCC lobes of peak negative 
(a) around the carbon atom are similar to those of 
the pure covalent bond. The kinetic energy density, as 
for C 3 H 4 , is concentrated in atom cores with little contri¬ 
bution from within the bond, and thus the bond’s polar 
character causes no observable change from that of the 
covalently bonded system. 

Next we consider SiH 4 , a nearly spherical molecule 
closely resembling a filled-shell atom in structure. It 
exhibits straightforward tetrahedral bonding, with sp^ 
hybridization of the silicon orbitals and covalent Si-H 
bonds. Fig. [^ shows a cut through a plane containing 
three of the atoms (H, Si, H) of the pseudo-molecule. A 
pair of hydrogen atoms is located above and below the 
plot plane as indicated by the open circles. 

An item of interest is the comparison between tks 
(c) and |VnJ^/n (a). Recalling that the von Weizsacker 
KED [Eq. (13)] is | Vn|^ / 8 n, we set the color scales of (a) 
and (c) to an exact 8:1 ratio so that a comparison of 
relative to t^w can be made. (For the other molecules, 
such a scheme wipes out almost all information about 
the gradient of the density, because tks is much larger 
than Tvw■) Here it is apparent that the Kohn-Sham KED 
reaches the von Weizsacker limit everywhere in the vicin¬ 
ity of a hydrogen atom. This seems reasonable in that 
each hydrogen atom has a single occupied orbital, and 
is in a sense a paradigm for the von Weizsacker limit in 
fermionic systems. The Laplacian for the Si-H bond (b) 
heavily emphasizes the H atom because of the more dis¬ 
persed nature of Si valence orbitals as compared to those 
of H. 

A final example from the test set is SiO, which features 
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FIG. 2. Functionals of the density for C 3 H 4 within the pseudopotential approximation, showing cut through the C-C-C bond 
axis and two hydrogens. In-plane atoms and bonds are shown as black disks and line segments; specific atoms are identified 
by labels. Out-of-plane ones shown as dashed lines and open disks. Contour levels for Laplacian (c) and pseudo-Laplacian 
of Eq. ( |15[ ) (d) are identical, with thick contour at zero. Contour levels for mCGA KED (f) are the same as those of the 
Kohn-Sham KED (e). 


a double bond that should be polar covalent given an 
electronegativity difference of 1.6. Fig. shows a surface 
plot of the Laplacian of the density for a cut through the 
pseudo-molecule Si-0 bond. Other quantities are avail¬ 
able for SiO in the supplementary material. The valence 
electrons that participate in this bond heavily favor oxy¬ 
gen, the more electronegative atom, leaving the silicon 
atom hypovalent.l^The zero contour (thicker contour) of 
is of interest for this system as it indicates a bond 
topology qualitatively different from the other cases. The 
orientation of the zero contour, crossing the bond perpen¬ 
dicularly to the bond axis, indicates that V^n > 0 in the 
bond center - specifically in the region between the two 


closed contours surrounding Si and O respectively. This 
indicates in QTAIM analysis that the bond is ionic; the 
maximum value of V^n = 0.194 is comparable to that of 
NaCl.l^ There is essentially a catastrophe in the topology 
of the zero contours whereby the Si-0 bond cannot be 
mapped to a polar-covalent geometry, such as the C-0 
bond in C 2 H 2 O 2 , without breaking and rejoining con¬ 
tours. 

Surface plots for the last two molecules of the AE6 test 
set (S 2 and C 4 H 8 ) are not shown in this paper but are 
available in the supplementary data.^ Although they re¬ 
peat themes already discussed for the other molecules, 
they have individual characteristics that should benefit 





































































(c) tks 


(d) TmGGA 


FIG. 3. Functionals of the density for pseudo-C 2 H 2 02, showing cut through the bond plane and an oxygen, carbon, and 
hydrogen atom. Details are the same as in Fig. 


from further investigation. The triplet ground state of S 2 , 
with a double bond and two lone pairs per atom, is simi¬ 
lar in structure to C 3 H 4 and C 2 H 2 O 2 . Nevertheless, the 
KED shows interesting structure near the valence shell 
peak and in the bond. Cyclobutane (C 4 H 8 ) is a cyclic 
molecule with a ring of four carbons and two hydrogens 
bonded to each. Unique to this system is the low-density 
region inside the carbon ring where the gradient of the 
density is zero but the Laplacian and KED are not. This 
topology is similar to that of the bond-center of a nearly 
dissociated molecule, and not found elsewhere (at equilib¬ 
rium geometry) in the test set. Such regions have been of 
interest for QTAIM analysi^^ and may provide a glimpse 
into how well approximated KED’s perform in predicting 
binding. 


C. Gradient expansion for the Laplacian 

The subfigure (d) of Eig. and (b) of Fig. show 
the pseudo-Laplacian V^n [Eq. (l^ ] which approximates 
the Laplacian in terms of the electron density, its gradi¬ 
ent and the kinetic energy density. Up to a small cor¬ 
rection proportional to |Vn| /n, this quantity is simply 
6 (Tifs — ttf); given that ttf is a power of the particle 


density, it interprets the Laplacian as roughly a measure 
of the difference between the kinetic energy and particle 
densities. As seen especially in Fig. our data sup¬ 
port this qualitative picture. The Laplacian (c) is posi¬ 
tive and large in the carbon pseudo-cores, precisely where 
the kinetic energy density (e) is largest and the charge 
density (a) is at a minimum; it is most negative in the 
bond regions where the situation is reversed. As a re¬ 
sult, V^n, plotted in (d) with the identical set of con¬ 
tours as V^n, captures the basic qualitative trends of 
this quantity, and on average, its relative magnitude in 
each bond. In contrast, the zero contour of V^n and 
V^n, shown as thicker black contours, have qualitatively 
different topologies. However, it seems reasonable to ex¬ 
pect that, given their qualitative similarity, they could 
produce similar results if used as parameters in a func¬ 
tional for an integrated quantity such as the exchange 
energy. Notably, the contour of V^n = 0 closely matches 
the shape of the 1/2-contour of the LOL, a close equiva¬ 
lent when density gradients are small . 1 ^ 

A check on the quality of this approximation can be 
obtained by the sum rule for V^n. Since it is an exact 
derivative, the integral of over the entire unit cell 
should be exactly zero. While the in tegra l for is 
zero to within round-off error, that for V^n ranges from 
about 0.1 hartree for SiH 4 to about 20 hartree for the 
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FIG. 4. Functionals of the density for pseudo-SiH 4 , showing cut through a plane containing a Si 
tetrahedrally bonded hydrogens. Details are the same as in Fig. 
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FIG. 5. Laplacian of the density for pseudo-SiO, showing cut through the bond axis. Details same as for Fig. j^c). 


largest molecules. This is a reflection of the very large 
difference between the integrated Kohn-Sham kinetic en¬ 
ergy and that of the Thomas-Fermi approximation. En¬ 
ergy densities can differ by several orders of magnitude 
in the pseudo-atom cores, an effect beyond the scale of 
our surface plots, but clearly shown in the log plots in 
Sec. HVEl 


D. The mGGA model for the kinetic energy density 

The final quantity for which we have made surface 
plots is the mGGA orbital-free KED.^^ It is shown for 
three molecules, subfigure (f) of Fig. and (d) of Fig. 
and[^ with contours and color scale that duplicate those 
of the Kohn-Sham KED. The agreement between the two 
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is generally not good. For C 3 H 4 (Fig. the mGGA, 
like the true KED, peaks in the pseudo-atom core, but is 
much smaller in magnitude. It is too large in the center of 
bonds, particularly the G^-C^ multiple bond. The most 
striking difference is the dramatic drop in magnitude in 
the mGGA in the region of peak valence charge con¬ 
centration surrounding each carbon atom. The shape of 
these zeroed-out regions correlates with the VSGG lobes 
of the Laplacian accentuating regions of peak density. 
The identical pattern shows up in C 2 H 2 O 2 (Fig. [^, with 
the KED zeroing out in VSGG regions for both oxygen 
and carbon, almost perfectly matching the contours of 
for the two lone oxygen pairs. This pattern occurs 
around the carbon atoms of C 4 H 8 , the two lone pairs 
of each S atom in S 2 and of the oxygen atom of SiO, 
indicating a global trend. 

SiH 4 , shown in Fig. is a case in which the mGGA 
works. In this case, much of the system is already very 
close to the von Weizsacker limit, which the mGGA is de¬ 
signed to capture exactly. Moreover, errors in the mGGA 
in different regions, such as in the Si atom core and near 
the hydrogen atom, almost exactly cancel, leading to a 
qualitatitively much better match of the mGGA to the 
exact KED than for the other five cases. (Notably, the 
Si pseudo-atom lacks the strong VSGG lobes associated 
with unusually low mGGA KED in G 3 H 4 and G 2 H 2 O 2 .) 


E. Plots through bond axes 


In this section, we focus on the quantitative compari¬ 
son of various models for the kinetic energy density, for 
which linear plots are convenient. We plot the enhance¬ 
ment factor Fs = t/ttf, which avoids excessive differ¬ 
ences in scale between atoms. We are also interested in 
the measure of electron localization a [Eq. (16)], that 
can also be thought of as the Pauli contribution to the 
enhancement factor. 

In Fig. we show Fs of several model KED’s for the 
pseudopotential approximation to the disulfur molecule 
S 2 as a function of displacement z from the molecule cen¬ 
ter along the bond axis. The focus is on a single sulfur 
pseudo-atom, marked by the black dot on the Fs = 0 
axis; the molecule has a mirror-symmetry plane through 
the bond center at z = 0. The Thomas-Fermi result is 
the horizontal line Fs = l. The von Weizsacker enhance¬ 
ment factor, TywlTTF = ^plS, is nearly zero in the bond 
region and again at the density peak associated with the 
lone pair behind the bond. The related expectation Uyw 
has an enhancement factor equal to lOg/S. It is negative 
in the covalent bond and the lone-pair behind the sul¬ 
fur atom, and positive in the pseudopotential core and 
asymptotically. 

Both the gradient expansion tcea and the more so¬ 
phisticated TraGGA are positive definite, in agreement 
with the required behavior of tks- Each is at a minimum 
in the bond and lone-pair regions, reach local maxima in 
the core and tend to 00 asymptotically. However, tks 
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FIG. 6. (color online) The enhancement factor Fs = t/ttf 
for the Kohn-Sham kinetic energy density (solid black line) 
and various orbital-free models, within the pseudopotential 
approximation, versus position along the bond axis 2 for S 2 . 
Also shown is Fs for u^w = n j The Thomas-Fermi result, 
trivially one, is shown as a dotted line. Location of the sulfnr 
atom noted by solid dot. 


is smooth and featureless, lacking the oscillatory struc¬ 
ture of the gradient and Laplacian of the density. The 
mGGA, where it differs from the GEA, does a slightly 
worse job in describing the Kohn-Sham value. In the 
lone-pair region around 2 = 3 a.u., it suffers from the 
extinction effect seen in the surface plots for G 3 H 4 and 
C 2 H 2 O 2 . Here, UywI ^tf < —1, equivalent to g < —0.3, 
which proves to be a significant criterion for this problem 
to occur in the mGGA. Overall, the mGGA fares better 
for S 2 than for other molecules, perhaps because this er¬ 
ror in its enhancement factor is cancelled by a reverse 
effect at the center of the double bond. The electron- 
localization measure a, not shown in Fig. is available 
in Fig. S-1 of the supplementary material.l^ 

Fig.0 shows enhancement factors for pseudo-SiO. As 
noted previously, this is the most polar molecule in the 
test set and gives a structural contrast to the more co¬ 
valent molecules. As such we focus on Tyw and Uyw as 
stand-ins for the gradient and Laplacian of the density, 
and related variables p and q, as compared to the Kohn- 
Sham KED. The gradient squared of the density {'^Tyw) 
does not vanish in the bond, as the density steadily in¬ 
creases from the Si valence shell to the O. The Laplacian 
Uyw) is positive at the center of bond, the QTAIM 
indication of ionic character. It is also instructive to plot 
the electron localization measure a, shown as the lighter 
(cyan) solid line in Fig. In the SiO bond, this measure 
approaches 0.5, equivalent to an ELF of 0.67, which is 
the value approached by the other double bonds in the 
test set. A more telling structure occurs behind the Si 
atom, where a falls nearly to zero over an extended re¬ 
gion of space. The value a ^ 0 (ELF of one) indicates that 
there is at most one occupied orbital so that t^s reaches 
von Weizsacker limit almost perfectly. It also coincides 


11 










z (a.u.) 


FIG. 7. (color online) The enhancement factor Fs=t/ttf 
for the Kohn-Sham and von Weizsacker KED’s, compared to 
the difference between the two, a, versus position along the 
bond axis z for pseudo-SiO. Also shown is Fs for u^w ~ V^n 
and the mGGArev discussed in Sec. |V G| dotted line indicates 
the Thomas-Fermi result. Location of each atom on bond axis 
noted by a solid dot. 



z (a.u.) 


FIG. 8 . (color online) The enhancement factor Fs = tIttf, 
plotted on a log scale for various kinetic energy densities ver¬ 
sus position along the carbon-carbon bond axis 2 for the C 3 H 4 
(propyne) pseudo-molecule. Location of atoms on axis noted 
by a solid dot. The Thomas-Fermi limit FIs = 1 is shown as 
dotted line. 


with an abnormally low minimum in tks- This probably 
is a reflection of the hypovalent character of Si in this 
molecule; however restricting the plot to the bond axis 
also eliminates the contribution of two 7r-bond orbitals 
to the KED. 

Fig.j^shows enhancement factors for the C 3 H 4 pseudo¬ 
molecule, for points through the axis joining the three C 
atoms and the on-axis terminal hydrogen (H^). We plot 
Fs on a log scale to focus on the situation at low densities, 
characterized by the carbon pseudo-atom cores and the 
asymptotic region far from the molecule. 

The asymptotic behavior of the Kohn-Sham and other 


KED’s is dominated by linear trend of log {Fs) to infin¬ 
ity far from the molecule {\z\ > 5). This is consistent 
with exponential decay of the charge density - and with 
ttf decaying more rapidly than any other model. 

The three orbital-free models shown - the von Weizsacker 
model, the GEA and the mGGA - have roughly the same 
decay constant, and for the most part match up quite 
well with the Kohn-Sham value. Interestingly, the GEA 
is the best predictor of tks, performing better than the 
mGGA almost everywhere. The von Weizsacker form al¬ 
most matches the Kohn-Sham case for the asymptotic 
edge near the lone hydrogen (2 > 5) - an indication that 
a single frontier orbital dominates the behavior of tks in 
this region. On the other edge of the bond axis (z<—5) 
Tks is roughly twice as large as t^w This area sees the 
intersection of three frontier orbitals, one from each of 
the three G-H bonds that form tetrahedrally off the cen¬ 
tral bond axis. This is enough to detach tks from the 
single occupied-orbital limit. 

An interesting story also occurs in the pseudopotential 
cores of the carbon atoms, with similar behavior seen 
for other atoms that have had core electrons replaced 
by pseudopotentials. Although this is arguably the least 
physical region of the molecule, it does represent one of 
rapidly varying low density, but negligible density gradi¬ 
ent, a topology that occurs in noncovalent bonds and the 
interstitial regions of solids. Here again the Kohn Sham 
KED is much larger than the Thomas Eermi value - as 
noted before, the charge and kinetic energy densities of 
our pseudopotential systems observe a kind of comple¬ 
mentarity, with one being large where the other is small. 
Of the three model KED’s, it is the GEA that repro¬ 
duces the KS value most accurately. The von Weizsacker 
model peaks at the edge of the core region where | Vn|^ is 
large, disappearing in the center of the pseudopotential 
core where it goes to zero. The GEA here closely follows 
which has a local maximum in the core and thus 
the correct qualitative behavior; surprisingly, the result 
is even quantitatively accurate. The mGGA trends more 
with TvWi and is severely deficient in magnitude. 

It is also useful for assessing approximate KED’s to 
plot the approximation to the electron-localization factor 
a obtained within a given model Tapprox'- 

^approx — {'^approx '^vw) i^TF ■ ( 1 ^) 

Focusing again on G 3 H 4 , which has the richest electronic 
structure of the test set, we plot in Fig. the a for sev¬ 
eral model KED’s on a log scale versus position z along 
the central bond-axis. For the exact Kohn-Sham t, we 
find three regions with a < 1 , an indication of electron lo¬ 
calization - the two carbon-carbon bonds and single ter¬ 
minating hydrogen atom. The other limit, a 1^1, occurs 
inside the pseudo-atom cores and asymptotically. The 
degree of localization of each small-a region is consistent 
with the character of that region. It is weakest (a ^0.5) 
for the triple bond between the second and third carbons, 
stronger (a ^ 0.3) for the single bond between the first 
two carbon atoms and extreme (a ^0.05) for the final hy- 
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FIG. 9. (color online) The electron localization measure 
ctapprox = {t — Tvw)/ttf, plotted ou a log scale for various ki¬ 
netic energy densities versus position along the carbon-carbon 
bond axis 2 : for the C 3 H 4 (propyne) pseudo-molecule. Loca¬ 
tion of atoms on axis noted by a solid dot. The Thomas-Fermi 
value a = 1 is shown as dotted line. The mGGArev is defined 
in Sec. lYl 


drogen atom where only a single orbital is occupied. As 
expected from the other figures shown, no approximate 
model does very well in these important situations. 

Two subtle order-of-limits issues come into play 
asymptotically. The approximate GEA and mGGA both 
do considerably worse in predicting the asymptotic trend 
of a in Fig.j^than they do the enhancement factor of tks, 
in Fig. a measures a difference between two models 
of r, and this difference is an order of magnitude smaller 
than the value of either model far from the molecule. It is 
thus a more sensitive probe of error in orbital-free mod¬ 
els. Secondly, while the GEA has the correct asymptotic 
behavior (although consistently three times too large), 
the mGGA has incorrect behavior as |z| —> 00 . To un¬ 
derstand this, note that a asymptotically can tend to 
anything from 0 to 00 . By the IP theorem, the numer¬ 
ator in Eq. (161 must vanish as the system tends to a 
one-electron state and tks^TvW, but the denominator 
also vanishes, as leaving the ratio undetermined. 

While the asymptotic value of a is one in the mGGA, for 
almost all the cases we have tested (SiO seems an excep¬ 
tion) the observed limit is infinity. This perhaps indicates 
only that ttf is an infinitely bad predictor of TpauU for 
a region in which the Thomas-Fermi approximation fails. 


V. ANALYSIS 

A. The GEA and asymptotic behavior of the KED 

It is worth analyzing in some depth what happens in 
the region of asymptotic decay far from the molecule, 
as demonstrated especially in Fig. and to some de¬ 
gree in Fig. It is striking that tks and tqea match 


each other almost perfectly in this limit, within 3% at 
higher densities and no more than 15% at the lowest 
densities we can obtain. Consequently V^n and its GEA- 
level approximation, V^n, also agree almost exactly for 
this region. This close agreement occurs for all systems 
studied, for example, for SiO as one either moves away 
from the hypovalent Si atom or from the nearly filled 
O atom. This is quite surprising since the GEA is de¬ 
signed for a completely different situation, that of the 
slowly varying electron gas, which is presumably unsuit¬ 
able for a classically forbidden region of space. Formally, 
the regime of validity of the slowly varying electron gas 
is for systems for which the inhomogeneity measures p 
and q are everywhere <?; 1. Obviously this criterion can¬ 
not be exactly met for a molecule, but one might expect 
that, in any extended region where these parameters are 
small, Tks should approach tgea- In fact, the opposite 
proves true: regions of space like covalent bonds, where p 
and q are consistently smallest, are where the GEA does 
the worst, while the classically forbidden asymptotic re¬ 
gion, where both p and q are much greater than one, is 
where it performs best. (To compare with the quanti¬ 
ties shown in Figs. and recall Tvw/ttf = 5p/3 and 
Uvw /ttf = ^0q/^-) Thus we have to conclude that some 
other phenomenon than the physics of the slowly-varying 
electron gas must explain the agreement asymptotically. 

It is not hard to find one, at least qualitatively. This re¬ 
gion is characterized by an exponential decay of the den¬ 
sity, n~exp (—2fcr), where k = gives the decay rate 
of the frontier orbitals, which have the highest eigenen- 
ergy, equal to the ionization potential I, and tunnel far¬ 
thest into the vacuum. As a result, the Kohn-Sham KED 
should behave as k^n, decaying at a rate proportional to 
the local particle density. In contrast the Thomas-Fermi 
KED varies as with kF . The enhancement 

factor Fg needed to correct t^f to the Kohn-Sham value 
then scales as k"^/kp, causing the exponential growth seen 
in Figs. 1^ and It is notable that the second-order gra¬ 
dient expansion reproduces this scaling behavior. The 
inhomogeneity variable p is equal to (k/kF)'^ for any ex¬ 
ponentially decaying particle density - and q is also, up to 
a correction due to curvature. The form of tgea [Eq. (|^] 
gives its enhancement factor the correct limiting behav¬ 
ior as r —>■ 00 . In contrast, the fourth-order correction 
has terms which scale like p^ = k^/kp, which blow up ex¬ 
ponentially as r—^oo. And a GGA, a closed expression 
summing over all orders of the gradient expansion, is not 
necessary to capture order-of-magnitude trends and can 
actually be less accurate than the second-order GEA. 

This is in stark contrast to what happens for the ex¬ 
change energy: the energy density associated with a sin¬ 
gle frontier orbital behaves asymptotically as (^)n while 
the LDA scales as kpn. Applying the second-order gra¬ 
dient expansion to the LDA creates an exchange energy 
density that scales incorrectly as {k'^/kp)n and a poten¬ 
tial that diverges exponentially. A GGA is needed to pro¬ 
duce an accurate exchange energy and a potential that 
is finite (if not with the correct 1/r form.) This contrast 
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between the correction needed for LDA exchange ^ 
and Thomas-Fermi KED ~ contradicts the conjoint¬ 
ness conjecture in its usual formulation - the same form 
of enhancement factor cannot be optimal for both cases. 


B. Revisiting the mGGA 

It is not hard to diagnose why the mGGA KED has 
difficulty modeling the molecular bond. The mGGA was 
tested primarily for closed-shell atoms and several model 
one-dimensional systems.!^ The most serious defects of 
the mGGA seen in the current study are associated with 
regions of of joint {p, q} space that these systems do not 
access. The issue of vanishing KED is strongly corre¬ 
lated with regions where p ~ 0 and q is negative and of 
the order of unity, a combination that does not happen 
with atoms.l^ A second problem is the mGGA’s large 
underestimate in the pseudopotential core where p ~ 0 
and q » 0; in atoms, a large g > 0 is associated with 
finite and normally large p, and occurs primarily in the 
asymptotic region far from the atom where both tend to 
OQ.HIl While both these errors lead to underestimates of 
the KE, the former is a failure to model the KED in the 
valence shell of an atom or molecule and should have a 
large impact on the model’s ability to predict molecular 
structure. 

The common thread here is the behavior of the KED as 
a function of q for values of p <C |q| and near to zero. In 
Fig. [^we plot the p = 0 limit of the KED enhancement 
factor, Fs{0,q), for several KED models as a function 
of the scale-invariant factor q. By definition, Fs = 1 for 
ttf, as shown by the dotted horizontal line. Likewise, 
the von Weizsacker KED for p = 0 is zero for all q. The 
Fs of the fourth-order gradient expansion approximation 
reduces to 1-I-(20/9)g-|-(8/81)q^. This is nearly indistin¬ 
guishable to the second-order gradient expansion, linear 
in q, because the fourth order coefficient is so small. The 
model of interest is the solid red line, that of Perdew and 
Gonstantin. It starts off with the gradient expansion and 
applies further constraints. First, the von Weizsacker 
bound requires that the enhancement factor be greater 
than Fyw, in effect greater than zero. For (7<0 the KED 
must transition fairly quickly from GEA-like behavior to 
zero, as the GEA breaks this constraint at g ^ —0.45. 
The second imposed limit is that the enhancement fac¬ 
tor goes to 1 + Fyw in the limit of large positive q, seen 
for example in our data in pseudo-atom cores, but not 
shown in Fig. [T^ 

The flaws in the mGGA seem to be caused by its imple¬ 
mentation of these constraints. The most important, the 
extinction of KED for negative q and small p, is clearly 
the result of a transition from tcea to that zeroes 
out the KED for q < — 0.30, a value of q achievable in 
the vicinity of an atomic lone-pair or a covalent bond, 
especially in a pseudopotential system. This transition 
scheme implicitly invokes an order of bounds as follows: 

tgea > tks > TvW ( 19 ) 



FIG. 10. (color online) Kinetic energy density enhance¬ 
ment factor Fs{p,q) for orbital free kinetic energy density 
models, plotted versus gradient-expansion parameter q for 
p — 0. Thomas-Fermi limit {Fs = 1) indicated by dotted 
line. Gray shaded region shows values of F disallowed by the 
von Weizsacker bound given by Fs = 0. 


That is, it seeks an interpolation between the two limiting 
cases, which leaves very little room for smoothing out the 
transition.!^ It makes more sense to try an interpolation 
above the two limiting cases, assuming a constraint 

Tks > ma.x{TGEA,Tyw), (20) 

demonstrated by the blue dashed curve in Fig. Such 
a transition is smoother and thus more physically appeal¬ 
ing, and has the effect of enhancing rather than reducing 
the KED in high-density, low-q, regions. A smoother 
transition to zero should also produce a smoother kinetic 
energy potential which is important for a self-consistent 
density functional minimization. Abrupt changes in the 
enhancement factor of a Laplacian-based density func¬ 
tional can be disastrous when taking functional deriva¬ 
tives with respect to V^n, since these involve derivatives 
of the density up to 

The mGGA model can also be improved by relaxing 
the large-g cutoff that it imposes. As seen especially in 
Fig.§ the mGGA clearly overcorrects for the regions 
where g 3> 1, in the pseudo-atom core and asymptoti¬ 
cally. In both cases, the second-order GEA is a better 
approximation to tks and there is less motivation for a 
GGA correction to it than is the case for exchange. 


C. Revision of the mGGA and application to atomization 
energy 


We propose to make a revised mGGA following two 
simple points: imposing the von Weizsacker lower bound 
by means of Eq. (20) and relying on the second-order gra¬ 
dient expansion otherwise. This satisfies the constraints 
for the two main limiting cases of the KED - that of delo¬ 
calized electrons with slowly-varying density and that of 
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strong electron localization, and otherwise keeps physi¬ 
cally reasonable behavior for regions of high inhomogene¬ 
ity - pseudo-atom cores and asymptotic decay. We first 
define a measure of electron localization 2 that depends 
upon the difference between GEA and vW enhancement 
factors for the KED: 

20 40 

2 : = Fgea - Fyw - 1 = - -^P ( 21 ) 

The factor z can be thought of as a poor-man’s electron 
localization factor - an orbital-free expression for the a 
used to describe electron localization in metaGGA’s and 
from which the ELF is constructed. We then look for the 
simplest possible asymptotic transition between Fqea 
and Fyw that imposes the von Weizsacker bound, which 
occurs for z < —1. Adapting a form recently used to 
construct a V^n-based exchange functioiP^ results in an 
enhancement factor 

F^GGArev = ^ ^ | ^ ^ _ ^(^)] | 

( 22 ) 

where H is the Heaviside step function. This is shown in 
Fig. [^as a function of q for p = 0. 

The limiting behavior of this correction can be charac¬ 
terized by three cases, roughly analogous to those defined 
by the ELF: 

1. If z —>■ 0, then both p and q must become small. 
The density is slowly varying, and close to the ho¬ 
mogeneous gas limit, typical of metallic bonds. In 
this case, Eg goes to the gradient expansion form: 

Fs ^ Fyw -I- 1 -I- 2 : = Fqea (23) 

2. If z < 0, this means that either q becomes nega¬ 
tive or p —>■ -|-oo with a finite q. In this case, Eg 
approaches the von Weizsacker limit: 

Fs ^ F,w + 0(l/z). (24) 

This is the proper description of a region with 
strong electron localization, such as a covalent 
bond. 

3. If z^O, we get the same result as for z small: 

Eg —>■ Fgea- (25) 

The primary situation for which this limit applies 
is an exponentially decaying density, for which 
q —>■ 00 and z —>■ 2 O 9 / 27 . 

The final case also describes a situation with q^l and 
finite p, seen here in pseudo-atom cores, and in the tran¬ 
sition between atomic shells in real atoms. 

Also of interest for real atoms is the limit q, z —>■ — 00 
which occurs near the nucleus and is caused by the 
cusp in the electron density. The functional derivative 
(5r/5n(r), used for the self-consistent determination of 


the charge density in OFDFT, must tend to Z/r near 
the nucleus so as to cancel the —Zjr contribution from 
the electron-nucleus potential. This behavior is exactly 
given by the functional derivative of and thus by 
TmGGArev as well. The leading correction to is of or¬ 
der I/g; its functional derivative is known to be finite, 
but it can cause a sizable error in the cusp of r at the 
nucleus. 


To evaluate the effects of this revised mGGA, we plot 
its enhancement factor for SiO in Fig. and its ap¬ 
proximation to a using Eq. (18 1 for propyne in Fig. 

As shown in the latter, the mGGArev by construction 
follows the GEA curve almost everywhere in space - 
except for regions of electron localization, where it en¬ 
hances the magnitude of the KED considerably over the 
GEA. It is thus an improvement over both GEA and 
mGGA. However the mGGArev overcorrects for situa¬ 
tions of strongest electron localization. For the single C^- 
G^ bond of propyne, with a small a of 0.3, the mGGArev 
gives a modest average overcorrection. It severely over¬ 
corrects for the most localized situations, where a < 0 . 1 : 
near the terminal H^ atom in propyne and behind the 
hypovalent Si atom in SiO. This problem may be ame¬ 
liorated by tinkering with the rate of transition between 
GEA and vW limits in Eq. (22) ~ in the current form 


(Fig. it is probably too slow. One region that shows 
little change from the mGGA is behind the G^ atom 
(z^—4) in Fig. This is not a region of electron local¬ 
ization since it feels the overlap of three neighboring G-H 
bond orbitals so the model has no criterion to correct for 
the error of the GEA. 


In Table HI we show errors with respect to the inte¬ 
grated Kohn-Sham KE averaged over the test set, as a 
measure of the overall quality of the models discussed in 
this paper. A net trend across all models is the underes¬ 
timation of the KE by roughly 10%. Unfortunately, by 
the virial theorem, the total KE is equal in magnitude 
to the total energy, which varies from 3.5 hartree for the 
valence shell of SiH 4 to 31 hartree for that of C 2 H 2 O 2 . 
Absolute errors in KE can thus be as large as several 
hartrees. While the second-order GEA is a modest im¬ 
provement over the Thomas-Fermi result, the mGGA, in 
attempting to address the limitations of the GEA, actu¬ 
ally loses some of the ground gained by it. The revised 
mGGA introduced here is more consistently an improve¬ 
ment. One situation in which it is not, SiH 4 , results in the 
maximum RE being three times the MARE, and an over¬ 
estimate, not an underestimate. As shown in Fig. this 
molecule is marked by a substantial region that is near 
the von Weizsacker limit, q:~0, for which the mGGArev 
overestimates the KED. This again indicates a need for 
further exploration of how to manage the transition from 
delocalized to localized electronic systems. 

To further characterize the quality of our revised 
mGGA, we calculate the atomization energies of the AE 6 
test set. This helps gauge the extent to which system¬ 
atic errors in the total energy are cancelled out in tak¬ 
ing energy differences. This is done not self-consistently. 
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TF 

GEA 

mGGA 

mGGArev 

VT84F 

MRE 

-0.162 

-0.112 

-0.124 

0.0021 

0.229 

MARE 

0.162 

0.112 

0.139 

0.0873 

0.229 

Max RE 

-0.202 

-0.159 

-0.178 

0.233 

0.550 


TABLE III. Mean relative error (MRE), absolute relative er¬ 
ror (MARE) and maximum relative error (Max RE) for vari¬ 
ous orbital-free estimates of the Kohn-Sham kinetic energy. 


System 

Exp. 

KS 

VT84F 

mGGA¬ 

rev 

mGGA 

GEA 

TF 

SiH4 

322.4 

315.9 

-178.1 

-14.9 

57.2 

-183.8 

-174.9 

S2 

101.7 

124.4 

140.9 

17.7 

-101.5 

-72.1 

-100.6 

SiO 

192.1 

205.4 

-4.8 

-4.5 

-169.3 

-97.6 

-213.7 

C 2 H 202 

633.4 

680.6 

476.6 

240.2 

-422.8 

-119.0 

-416.6 

G 3 H 4 

704.8 

726.6 

581.9 

572.3 

24.2 

115.9 

35.7 

G 4 H 8 

1149 

1175.3 

1072.4 

811.8 

142.6 

96.0 

-53.3 

MAE 

- 

23.0 

182.1 

246.8 

595.5 

560.7 

671.1 


TABLE IV. Atomization energies for the AE6 test set in 
kcal/mol. Shown are experimental values from Ref. 1601 self- 
consistent Kohn-Sham results, and results of orbital-free mod¬ 
els evaluated with the Kohn-Sham density. Also shown is the 
mean absolute error with respect to experiment. 


using conventional PBE Kohn-Sham densities and bond 
lengths (Table|l]). The results are shown in Table IV First 
we note how far the Thomas-Fermi atomization energy is 
from experiment, with an MAF ten times worse than the 
LDA and over thirty times worse than the PBF Kohn- 
Sham models. It almost always fails to predict binding, 
at best giving a marginal binding energy. The second- 
order GEA does provide a modest improvement over the 
TF case, but again shows severe under-binding. By re¬ 
specting the von Weizsacker lower bound, the mGGA 
ought to significantly improve GEA atomization ener¬ 
gies. Instead it performs worse for the majority of the 
test set, and in some cases worse than Thomas-Fermi. 
In constrast, the mGGArev does show the expected im¬ 
provement over Thomas-Fermi and GEA. It binds all but 
one molecule, SiH 4 , the standout worst case in Table [ITTl 
and the one case that the mGGA binds. On average it 
removes 60% of the AE error of the TF and for one or two 
systems almost approaches the LDA in quality. However 
its MAE is still an order of magnitude worse than the 
PBE and a factor of three worse than the LDA. 

To put these results in perspective, we perform cal¬ 
culations for the VT84F,E^ a nonempirical GGA for the 
kinetic energy. This applies the key constraints of the 
mGGA - respecting the gradient expansion in the small- 
p limit and requiring the von Weizsacker constraint for 
all p] in addition it enforces the non-negativitji^ of the 
Pauli potential, STpauii/Sn{r) > 0. The VT84F total ki¬ 
netic energy (Table III) is by a large margin the least 
accurate of all models considered, including the Thomas- 
Fermi model. However it has the overall best prediction 


of atomization energies (Table [rv| ), and fails significantly 
only for SiH 4 . It may be hard to enforce both the GEA 
and the constraint r > with only access to |Vn|^ as 
a variable and not overestimate the total KE. However 
enforcing constraints on the potential - an infinitesimal 
energy difference - seems to help for predicting accurate 
finite-energy differences. It is reassuring that the sim¬ 
ple metaGGA we present here is comparable in quality 
to the VT84F without (as yet) taking the potential into 
consideration. 


VI. DISCUSSION AND CONCLUSIONS 

We present highly converged DFT calculations for the 
AEG test set of molecules, within a plane-wave pseudopo¬ 
tential approach. We use these to visualize the Kohn- 
Sham kinetic energy density and related quantities that 
are ingredients of modern DFT’s, specifically metaGGA 
models for the exchange-correlation energy, and orbital- 
free models for the KED. By providing a highly accurate 
map between density and kinetic energy density for phys¬ 
ically reasonable model systems, our data enables the use 
of visualization techniques employed in the qualitative 
analysis of electronic structure to test approximations 
to this critical ingredient for DFT. The pseudopotential 
method works especially well in characterizing the classi¬ 
cally forbidden region far from nuclei, and is reasonable 
in its description of bonds; its main limitation is the loss 
of knowledge of the core region, most importantly, the 
character of the core-valence transition that plays a key 
role in determining bond lengths. 

The choice of the AEG test set does not break new 
ground in visualization of electronic structure, but does 
an excellent job of illustrating many of the lessons learned 
from QTAIM and other visualization approaches, par- 
ticurly the role of in understanding valence elec¬ 
tronic structure and the KED in measuring electron lo¬ 
calization. The SiO molecule is perhaps of most interest 
structurally, given the relationship between the hypova- 
lent character of Si and the strong indication of electron 
localization in the Si valence shell; also of interest is the 
identihcation of the bond as ionic rather than polar co¬ 
valent by QTAIM criteria. 

A major finding of the paper is the surprising success 
of the gradient expansion expression for the Kohn-Sham 
KED. The gradient expansion approximation V^n ^ 
6{tks~'''tf) used in modern metaGGA’s is at least qual¬ 
itatively very good - V^n to some degree picks up the 
complementary behavior of the kinetic and particle den¬ 
sities, and detects regions where one is larger than the 
other. Rather surprisingly, this approximation is the 
most accurate in the lowest density regions, in the clas¬ 
sically forbidden regions far from nuclei. This is because 
it has the exact asymptotic behavior with respect to dis¬ 
tance from nuclei and not too bad quantitative values for 
all systems considered. 

The asymptotic exactness of the GEA, although not 
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news, is worthy of note since it points out the limitations 
of the idea of conjointness between exchange-correlation 
and kinetic energy functionals, both as a conjecture and 
as a design philosophy. Lessons learned in designing func¬ 
tionals for the former case do not necessarily transfer over 
to the latter. The very different behavior of the gradi¬ 
ent expansion for exchange in the asymptotic limit ne¬ 
cessitates a fundamentally different functional form for 
exchange energy GGA’s and kinetic energy GGA’s. The 
gradient expansion of the former must be controlled by 
some form of cutoff at large values of reduced density 
gradient p while that of the latter, as best we can see, is 
better off mostly untouched. 

A second point underscores the difficulty in building 
orbital-free models of the KE - the gradient expansion 
behaves worst in describing “slowly varying” regions of 
space - where the inhomogeneity parameters p and q used 
to describe it are small. For the KE density, there seems 
not to be a good “semilocal” approximation for real sys¬ 
tems - one cannot rely on p and q being small locally 
to predict that the gradient expansion should hold lo¬ 
cally, in contrast with the XC energy density. When 
taken separately, exchange and correlation energy densi¬ 
ties have similar problems to those we see here for the 
KED; however, there is a notable cancellation of error 
between the two that makes semilocal approximations 
work better than might be expected.!^ What the KED 
lacks then is a companion mechanism such as correlation 
by which deviations from the GEA can be cancelled out. 
This failure does not contradict the idea of the gradi¬ 
ent expansion. The limit in which it is exact is that of 
globally small p and q, with the result of delocalized elec¬ 
tronic orbitals almost everywhere, a condition that is not 
met by any molecule. 

The other major finding of this paper relates to the 
Perdew-Constantin mGGA model of the kinetic energy 
density. We have found a number of problems which de¬ 
grade its performance with respect to the Thomas-Fermi 
model. Its description of the KED in regions of high in¬ 
homogeneity and low density are less effective than the 
simpler second-order GEA. More importantly, it is sub¬ 
ject to an “extinction” effect for large negative values of 
the reduced Laplacian q, causing it to plummet to zero in 
regions of covalent bonding. This effect is caused by the 
particular form used in the imposition of the constraint 
which becomes important in regions of electron 
localization, such as covalent bonds. It is aggravated by 
the use of pseudopotentials, which exaggerate the magni¬ 
tude of negative-V^n or VSGG regions in comparison to 
their all-electron counterparts. The correlation between 
this effect and electron localization seems responsible for 
the poor binding seen with this model. The formation of 
bonds can reduce electron localization and thus reduce 
the extinction effect relative to the isolated atom case, 
leading to a lack of error cancellation in taking energy 
differences. Notably, the cases in which the mGGA gives 
improved binding energies, SiH 4 and C 4 H 8 , are the ones 
with exclusively single bonds and thus roughly the same 


degree of electron localization in molecule and atom. 

This work points to several avenues of future research. 
The mGGArev form we propose for the KED is the sim¬ 
plest, not best, form that can ht the constraints imposed 
in Sec. |V] and should perhaps be used not as a finished 
functional but as an indication of how to proceed in devel¬ 
oping one. Particularly, the use of an “orbital-free ELF”, 
using derivatives of the density to approximate the ELF 
and its ability to distinguish between different kinds of 
bonds, seems worthy of further investigation. However, 
in its current form, our model regresses on the mGGA’s 
capacity to handle covalently bonded hydrogen atoms 
and other situations of nearly perfectly localized elec¬ 
trons. Notably, our proposed constraint, that t>tgea 
when T^Tvw, is not universal - it fails for the Is shell of 
atoms, as shown in Ref. [23l Not surprisingly, the mGGA 
functional, with r < tgea in this limit was arrived at 
partly through the consideration of this case. However 
our constraint does appear to be valid for any other shell 
of an all-electron atom - and it is responsible for our 
current revision’s relative success in predicting binding 
energies of the AEG test set. Any more sophisticated 
model will thus have to ameliorate somehow the prob¬ 
lems for hydrogen encountered by the mGGArev while 
keeping its nice features for bonding. 

A second notable issue is the large deviation of V^n 
and tks obtained with pseudopotentials from their all¬ 
electron values just inside the pseudopotential cutoff ra¬ 
dius. As noted earlier, the resulting exaggeration of the 
negative value of in VSCC’s contributes to the fail¬ 
ure of the mGGA to predict binding in pseudopotential 
systems. But, given the sudden switching behavior that 
the mGGA shows for negative q (Fig. [^, it is quite 
possible that with the smaller q values of all-electron 
systems, this effect would not be an issue. This raises 
the question of how much pseudopotentials that have 
been constructed to match valence density alone can be 
trusted in metaGGA or OFDFT applications that rely 
upon variables that are more sensitive to changes in elec¬ 
tronic structure. 

In the big picture, the ability to get orbital-free DFT’s 
that are competitive with Kohn-Sham approaches re¬ 
mains a challenge. However, some progress towards 
functionals useful in extreme situations where Kohn- 
Sham approaches are impractical may yet be done with 
metaGGA’s working with semilocal properties of the den¬ 
sity. Visualization can be a useful tool in this process, 
fruitfully bringing together strands of qualitative and 
quantitative thinking about electronic structure. 
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